{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Extraction files"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## ICESat 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import glob\n",
    "import os\n",
    "import pandas as pd\n",
    "import geopandas as gpd\n",
    "from pathlib import Path\n",
    "import h5py\n",
    "import re\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.dates as mdates\n",
    "import time\n",
    "from astropy.time import Time\n",
    "import datetime\n",
    "import cartopy.crs as ccrs\n",
    "import cartopy.io.img_tiles as cimgt\n",
    "import statistics\n",
    "%matplotlib widget\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##### load files\n",
    "## set the directory\n",
    "data_home = Path('/home/jovyan/ICESat_water_level/extraction/download/')\n",
    "## list them up and check them\n",
    "files= list(data_home.glob('*.h5'))\n",
    "### bounds of Tonle Sap Lake\n",
    "sp_ex = [103.643, 104.667, 12.375, 13.287]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Raw files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### I made a function for stacking the information\n",
    "def alt13_to_df_beams(filename):    \n",
    "    f = h5py.File(filename, 'r')\n",
    "    rgt = str(filename).split(\"_\")[5][0:4]\n",
    "    cycle = int(str(filename).split(\"_\")[5][4:6])\n",
    "    version = int(str(filename).split(\"_\")[5][6:8])\n",
    "    ymd = str(filename).split(\"_\")[4][0:8]\n",
    "    ymd_trans = datetime.datetime(int(ymd[0:4]),int(ymd[4:6]),int(ymd[6:8]))\n",
    "    date = ymd_trans.strftime(\"%Y-%m-%d\")\n",
    "    beam_lst = list(f)[2:-1]\n",
    "    alt13_df = pd.DataFrame()\n",
    "    for beam in beam_lst:\n",
    "        f_beam = f[beam]\n",
    "        lat = f_beam['segment_lat'][:]\n",
    "        long = f_beam['segment_lon'][:]\n",
    "        ws = f_beam['ht_water_surf'][:]\n",
    "        ws_sd = f_beam['stdev_water_surf'][:]\n",
    "        ws_err = f_beam['err_ht_water_surf'][:]\n",
    "        ortho = f_beam['ht_ortho'][:]\n",
    "        wd = f_beam['water_depth'][:]\n",
    "        df_beam = pd.DataFrame({'Beam': beam ,'RGT':rgt,'Cycle': cycle, 'Date':date, 'Date_num':int(ymd),'Ver.':version,\n",
    "                                 'Latitude':lat,'Longitude':long,'SurfaceH':ws,\n",
    "                                'SH_SD':ws_sd, 'SH_error':ws_err,'OrthoH':ortho,\n",
    "                                'WaterD':wd})\n",
    "        alt13_df = alt13_df.append(df_beam, ignore_index = True)\n",
    "        \n",
    "    return alt13_df\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a_tracks = pd.DataFrame()\n",
    "for ff in files:\n",
    "    a_tracks = a_tracks.append(alt13_to_df_beams(ff), ignore_index = True)\n",
    "a_tracks.sort_values(by=['Date_num'], inplace=True)\n",
    "print('done')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "icesat2 = a_tracks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "icesat2.to_csv('/home/jovyan/ICESat_water_level/extraction/data/icesat2.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Without outliers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### outliers\n",
    "def out_rmv(df,var):\n",
    "    ### lower (25%)\n",
    "    lwq = df[var].quantile(0.25)\n",
    "    ### upper (75%)\n",
    "    hwq = df[var].quantile(0.75)\n",
    "\n",
    "    ### IQR\n",
    "    iqr = hwq - lwq\n",
    "\n",
    "    lw_out = lwq-1.5*iqr\n",
    "    hw_out = hwq+1.5*iqr\n",
    "\n",
    "    ###LSHoutlier > q0.75 + 1.5 ∗ IQR or LSHoutlier < q0.25 − 1.5 ∗ IQR\n",
    "    return df.loc[(df[var] >= lw_out) & (df[var] <= hw_out)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### I made a function for stacking the information\n",
    "def alt13_a_out(db,var2):\n",
    "    alt13_df = pd.DataFrame()\n",
    "    date_list = db['Date_num'].unique()\n",
    "    for prd in date_list:\n",
    "        s_db = db.loc[db['Date_num']==prd]\n",
    "        beam_lst = db['Beam'].unique()\n",
    "        for beam_type in beam_lst:\n",
    "            df_select = s_db.loc[s_db['Beam']==beam_type]\n",
    "            alt13_df = alt13_df.append(out_rmv(df_select,var2), ignore_index = True)\n",
    "    return alt13_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "icesat2_out=alt13_a_out(icesat2,'SurfaceH')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "icesat2_out.to_csv('/home/jovyan/ICESat_water_level/extraction/data/icesat2_out.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Average"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### I made a function for averaginf the height\n",
    "def alt13_mean(db,var2):\n",
    "    mean_df_list = pd.DataFrame()\n",
    "    date_list = db['Date_num'].unique()\n",
    "    for prd in date_list:\n",
    "        s_db = db.loc[db['Date_num']==prd]\n",
    "        av_height = statistics.mean(s_db[var2])\n",
    "        mean_df = pd.DataFrame({'Av_level':av_height,'Date':s_db['Date'].unique(), 'Date_num':s_db['Date_num'].unique()})\n",
    "        mean_df_list = mean_df_list.append(mean_df, ignore_index = True)\n",
    "    return mean_df_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "icesat2_av = alt13_mean(icesat2_out,'SurfaceH')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "icesat2_av.to_csv('/home/jovyan/ICESat_water_level/extraction/data/icesat2_av.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## ICESat 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##### load files\n",
    "## set the directory\n",
    "data_home = Path('/home/jovyan/ICESat_water_level/extraction/icesat/')\n",
    "## list them up and check them\n",
    "files= list(data_home.glob('*.H5'))\n",
    "### Spatial Bounds: \n",
    "tsl_sp_ex = [103.643, 104.667, 12.375, 13.287]\n",
    "### load the water mask\n",
    "tsl_wm = gpd.read_file('/home/jovyan/ICESat_water_level/extraction/shp/tsl_sample_dis.shp') ## water mask: Tonle Sap Lake"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def glah14_to_df(filename,bounds):\n",
    "    ## sp_ex = [103.643, 104.667, 12.375, 13.287]\n",
    "    ## Bounds are [Longitude_West, Longitude_East, Latitude_South, Latitude_North]\n",
    "    f = h5py.File(filename, 'r')\n",
    "    lat = f['Data_40HZ']['Geolocation']['d_lat'][:]\n",
    "    lon = f['Data_40HZ']['Geolocation']['d_lon'][:]\n",
    "    elev = f['Data_40HZ']['Elevation_Surfaces']['d_elev'][:]\n",
    "    sec = f['Data_40HZ']['Elevation_Corrections']['d_satElevCorr'][:]\n",
    "    scf = f['Data_40HZ']['Quality']['sat_corr_flg'][:]\n",
    "    satndx = f['Data_40HZ']['Quality']['i_satNdx'][:]\n",
    "    dem = f['Data_40HZ']['Geophysical']['d_DEM_elv'][:]\n",
    "    date_gps = f['Data_40HZ']['DS_UTCTime_40'][1] + 630763213 \n",
    "    date_ymd = Time(date_gps,format='gps').datetime.strftime(\"%Y-%m-%d\")\n",
    "    \n",
    "    glah14_df = pd.DataFrame({'Latitude':lat,'Longitude':lon,'Elevation':elev,\n",
    "                            's_El_Corr':sec, 's_Corr_f':scf,'in_sat':satndx,\n",
    "                            'DEM':dem,'Date':date_ymd})\n",
    "    #### Subsetting\n",
    "    glah14_df_subset = glah14_df.loc[(glah14_df['Longitude']>=bounds[0]) \n",
    "                          & (glah14_df['Longitude']<=bounds[1])\n",
    "                          & (glah14_df['Latitude']>=bounds[2])\n",
    "                          & (glah14_df['Latitude']<=bounds[3])]\n",
    "    return glah14_df_subset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test1=glah14_to_df(files[0],tsl_sp_ex)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test2=glah14_to_df(files[9],tsl_sp_ex)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test1.append(test2)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "is1_pd = pd.DataFrame()\n",
    "for kk in range(0,5):\n",
    "    is1_pd = is1_pd.append(glah14_to_df(files[kk],tsl_sp_ex), ignore_index = True)\n",
    "    print(kk)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(is1_pd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#### convert to geoDF and clip it. \n",
    "is1_pd = pd.DataFrame()\n",
    "for num in range(1,len(files)):\n",
    "    gdf = gpd.GeoDataFrame(glah14_to_df(files[num],tsl_sp_ex))\n",
    "    gdf.set_geometry(\n",
    "        gpd.points_from_xy(gdf['Longitude'], gdf['Latitude']),\n",
    "        inplace=True, crs='EPSG:4326')\n",
    "    gdf_wm = gpd.clip(gdf, tsl_wm)\n",
    "    is1_pd = is1_pd.append(gdf_wm)\n",
    "    print(num)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [],
   "source": [
    "icesat1 = is1_pd\n",
    "\n",
    "icesat1.sort_values(by=['Date'], inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [],
   "source": [
    "icesat1.to_csv('/home/jovyan/ICESat_water_level/extraction/data/icesat1.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### remove outliers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [],
   "source": [
    "### I made a function for stacking the information\n",
    "def gl14_a_out(db,var2):\n",
    "    gl14_df = pd.DataFrame()\n",
    "    date_list = db['Date'].unique()\n",
    "    for prd in date_list:\n",
    "        s_db = db.loc[db['Date']==prd]\n",
    "        gl14_df = gl14_df.append(out_rmv(s_db,var2), ignore_index = True)\n",
    "    return gl14_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [],
   "source": [
    "icesat1_out = gl14_a_out(icesat1,'Elevation')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [],
   "source": [
    "icesat1_out.to_csv('/home/jovyan/ICESat_water_level/extraction/data/icesat1_out.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Average values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [],
   "source": [
    "### I made a function for averaginf the height\n",
    "def gl14_mean(db,var2):\n",
    "    mean_df_list = pd.DataFrame()\n",
    "    date_list = db['Date'].unique()\n",
    "    for prd in date_list:\n",
    "        s_db = db.loc[db['Date']==prd]\n",
    "        av_height = statistics.mean(s_db[var2])\n",
    "        mean_df = pd.DataFrame({'Av_level':av_height,'Date':s_db['Date'].unique()})\n",
    "        mean_df_list = mean_df_list.append(mean_df, ignore_index = True)\n",
    "    return mean_df_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [],
   "source": [
    "icesat1_av = gl14_mean(icesat1_out,'Elevation')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [],
   "source": [
    "icesat1_av.to_csv('/home/jovyan/ICESat_water_level/extraction/data/icesat1_av.csv')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
